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Abstract 

Tlie aim of tlie researcli article is to develop a representative volume element using finite elements to study 
the buckling stability of graphene/polymer nanocomposites. Research work exploring the full potential of 
graphene as filler for nanocomposites is limited in part due to the complex processes associated with the 
mixing of graphene in polymer. To overcome some of these issues, a multiscale modeling technique has been 
proposed in this numerical work. Graphene was herein modeled in the atomistic scale, whereas the polymer 
deformation was analyzed as a continuum. Separate representative volume element models were developed for 
investigating buckling in neat polymer and graphene/polymer nanocomposites. Significant improvements in 
buckling strength were observed under applied compressive loading when compared with the buckling stability 
of neat polymer. 
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Background 

In the recent past, graphene has emerged as a potential 
candidate for developing nanocomposites with improved 
properties [1,2]. The experimental characterization of gra- 
phene/polymer nanocomposites is a challenging process, 
and hence, computational approaches for predicting the 
behavior of such materials have also extensively been 
employed. Various multiscale models are available in the 
literature for predicting the properties of carbon nanotube 
(CNT) -based nanocomposites [3-5], but very few models 
have been presented to study graphene nanocomposites. 
For example, Cho et al. [6] developed a numerical model 
in conjunction with Mori-Tanaka approach to study the 
elastic constants of randomly distributed graphene in 
polymer. Awasthi and his team [7] investigated the load 
transfer mechanism between polyethylene and graphene 
sheets. Montazeri and Tabar [8] developed a finite element 
(FE)-based multiscale model to investigate the elastic con- 
stants of graphene-based nanocomposites. 
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Buckling in isolated graphene sheets was modeled by 
several researchers [9-11]. However, budding stability 
of graphene/polymer nanocomposites was only 
reported by Rafiee et al. [12]. Using an experimental 
and analytical approach, up to 50% and 32% improve- 
ment in the buckling stability of nanocomposites was 
reported respectively. In the analytical approach, an 
Euler buckling formulation was employed, and elastic 
properties required in the Euler equation were esti- 
mated by experimental means. The discrepancies be- 
tween the two buckling stabilities were attributed to 
scaUng issues. 

It is well established that the reinforcement of polymer 
with graphene increases the elastic modulus of the material 
which further improves buckling stability. The aim of this 
study is to propose a numerical model which can estimate 
the increase in buckling stability with different volume 
fractions of graphene and can further be extended to com- 
plex shapes and structures. 

It has been reported that achieving a uniform disper- 
sion of two-dimensional graphene sheets in polymer is 
more challenging compared to the mixing of one- 
dimensional CNT. Moreover, the application of nano- 
composites is not limited to simple structures, and the 
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comprehension of material behavior in complex struc- 
tures is restricted when employing experimental and 
analytical methods. Consequently, research efforts are 
increasingly focused on numerical approaches. To 
overcome some of the limitations that exist in experi- 
mental and analytical work, a multiscale representative 
volume element (RVE) is proposed in this paper to in- 
vestigate buckling phenomena in graphene/polymer 
nanocomposites under the assumption that graphene is 
uniformly distributed in the polymer. To the knowledge 
of the present authors, no numerical model has been 
reported yet to study the effect of graphene on the 
budding strength of nanocomposites. In the proposed 
technique, graphene was modeled in the atomistic 
scale, whereas polymer deformation was analyzed as a 
continuum. 

Methods 

Finite element modeling of RVE 

In this paper a finite element technique was employed 
in conjunction with molecular and continuum 
mechanics to simulate budding in graphene/polymer 
nanocomposites. In the proposed RVE, the polymer, 
epoxy in this case, was modeled as a continuum phase 
whereas the deformation in graphene was evaluated 
using an atomistic description. Nonbonded interac- 
tions were considered as the load transfer mechanism 
or interphase between polymer and graphene. FE 
modeling was performed in the ANSYS (Version 13) 
software environment (ANSYS Inc., Canonsburg, PA, 
USA). 

Atomistic model for graphene 

To model graphene in the proposed RVE, it was 
assumed that graphene behaves like a space frame struc- 
ture in which the covalent bonding between C-C atoms 
constitutes the load-bearing element while atoms act like 
a joint. A molecular-mechanics-based approach was 
employed to estimate properties for graphene in the 
atomistic scale. In the molecular form, potential energy 
of graphene can be represented as 



energies of a beam element as described by Equations 2, 
3, and 4. 
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In Equation 1 Vbonded represents the bonded inter- 
action between atoms and is considered primarily re- 
sponsible for maintaining the structural integrity. The 
potential energy term ' Vnonbonded' includes electrostatic 
forces, e.g., van der Waals forces, which are considered 
weak in nature and can generally be neglected in com- 
parison to bonded interactions. Li and Chou [13] 
established a correlation between interatomic molecu- 
lar potential energies and corresponding strain 
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where /C (938 kcalmol'^A"^), Kq (126 kcalmol'^rad"^), 
and (40 kcalmoP^rad"^) are the bond stretching, 
bond bending, and bond torsional resistance force con- 
stants. E represents the Young's modulus; A, the cross 
sectional area; /, the moment of inertia; G, the shear 
modulus; /, the polar moment of inertia of the beam 
element, respectively. BEAM4 elements were used in 
the ANSYS software environment to model graphene 
by connecting nodes, and material properties for those 
elements were estimated with the help of Equations 2, 
3, and 4. 

Continuum model for polymer 

The volume fraction of graphene in polymer ranges 
commonly up to 10%. Most of the material volume is 
therefore occupied by polymer. Simulating the polymer 
phase on the atomistic scale would require large efforts 
in dealing with large numbers of degrees of freedom as 
well as substantial computational cost. Therefore, as a 
reasonable compromise, the polymer phase was modeled 
as a continuum, and three-dimensional SOLID45 ele- 
ments were used for meshing the geometry. Epoxy with 
a Young's modulus of 3.4 GPa and a Poissons ratio of 
0.42 was considered as the polymer material in the 
present work. 

Interphase between graphene and polymer 

In this research paper, nonbonded interactions were 
considered as load transfer mechanisms between the 
polymer and graphene. Naturally, load transfer between 
nonfunctionalized graphene sheets and polymer takes 
place through van der Waals interactions. The Lennard 
Jones '6-12' potential given in Equation 5 was employed 
to estimate the properties for the interface region, which 
is depicted in the schematic in Figure 1. 
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where r is the atomistic distance between atoms, 0 (0.34 
nm) is the hard sphere radius, and y (0.556 kcal/mole) is 
the potential well depth. In the finite element environ- 
ment, these nonbonding interactions were modeled with 
the help of LINKS truss elements. As shown in 
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Equation 6, properties for those truss elements were 
estimated by comparing the classical continuum strain 
energy of truss elements with the Lennard Jones '6-12' 
potential Note that the truss model described herein 
was earlier employed by Li and Chou [14] to simulate 
the interaction between CNT and epoxy polymer. 
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where 7?Eq is the initial undeformed length, r is the 
deformed length, A is the cross sectional area, and E{r) 
is the Young's modulus of the truss element. 

Eigenvalue buckling analysis 

A linear analysis for mode one buckling was performed 
in this work, which is associated with the computation 
of a bifurcation load and corresponding buckling mode. 
The analysis in the finite element environment was 
divided into two sections, i.e., a pre-buckling and a post- 
buckling analysis [15]. 

Pre-buckling analysis 

The pre-buckling analysis was performed in the finite 
element environment to compute the reference stres- 
ses S"' within each element. These computed reference 
values are then used for estimating the geometric 
stiffness matrix Kg, which is needed in the post- 
buckling analysis. To initiate the analysis in the finite 
element model, the general matrix form equation pro- 
vided in Equation 7 was employed to estimate nodal 
displacements [Q] when the structure is subjected to 
the reference unit load Pj^, In Equation 7 the total 
stiffness matrix [Kq] is the assembled global stiffness 
matrix for the multiscale structure composed of 
beam, truss, and solid elements. Conventional 



assembly techniques can be used to obtain the global 
stiffness matrix. 



[Ko][Q] = [Pr] 



(7) 



Equation 8 represents the general solution that is 
obtained in the finite element analysis after solving the 
assembled matrices defined in Equation 7. 



{Q^} = [Kor'[PR] 



(8) 



The nodal displacement vectors {Q''} obtained from 
Equation 8 are further post-processed to compute strains 
and corresponding stress values 5'* within each element 
using Equation 9, in which matrix [D] represents the elas- 
ticity matrix (material property) while matrix [B] is a func- 
tion of approximating polynomial or shape functions. 



{S*} = [D][B]{Q*} 



(9) 



Post-buckling analysis 

The post-buckling analysis corresponds to a general 
eigenvalue problem as defined by Equation 10, in which 
A and v are the load factor and the eigenvector of displa- 
cements. The geometric stiffness matrix used in 
Equation 10 is defined with the help of Equation 11, 
where F is composed of derivatives of shape functions 
and 5' is a function of the stress 5* estimated in the pre- 
buckling analysis. 
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By solving Equation 10, the lowest possible value for 
the load factor A^in can be determined, which is then 
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23.79 nm 



Figure 2 Schematic of model and boundary conditions for estimating Young's moduli of developed RVE. 



employed to estimate the buckling force given by Equa- 
tion 12. 

Bucklingforce = A^/„ x Reference load (12) 
Results and discussion 

The proposed RVE was employed to understand the 
buckling behavior of graphene in a polymer matrix when 
graphene was assumed to be uniformly distributed. A 
second separate RVE structure with equal dimensions 
was developed using finite element modeling to study 
the buckling in neat (homogenous) polymer. 



To validate the results of the proposed RVE, an alter- 
native method (indirect method) was also employed in 
this paper. The proposed RVE model with boundary 
conditions as shown in Figure 2 was first used to esti- 
mate the Young's modulus of the RVE as a whole for 
various volume fractions of graphene. The estimated 
material properties (£) were then used in defining the 
homogenous rectangular plate with dimensions of the 
proposed RVE model. The buckling loads were calcu- 
lated for this homogenous rectangular plate shown in 
Figure 3a and were compared with the direct approach 
defined in Figure 3b. 



Compressive force 



Compressive force 



E 
c 

CO 

eg 



Neat 
polymer 



12 nm 



E 
c 

1^ 

CO 
CVJ 



Polymer 

I 

12 nm 



Graphene 
sheet 



WWW \ W\W" 

Constraint Constraint 
(a) (b) 

Figure 3 Boundary conditions and dimensions for (a) neat polymer model and (b) multiscale graphene/polymer model. 
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The analytical formulation to estimate the buckling load 
for a rectangular plate is given in Equation 13, where is 
the applied unidirectional compressive load, w is the dis- 
placement in the outward normal direction, and D is flex- 
ural rigidity (i.e., a function of the Young's modulus). 
Boundary conditions were kept identical in all RVE mod- 
els. Only the value of flexure rigidity was varied, which 
afforded the proposed RVE model the capability of yield- 
ing E for different graphene volume fractions. 

DvM^.y)^Pxg^ = o (13) 

The boundary conditions along with dimensions for 
the proposed RVE model are shown in Figure 3. The 
thickness of graphene in the atomistic scale and for 
epoxy as a continuum phase was kept constant at 0.344 
nm, whereas the thickness of the interphase was kept at 
0.172 nm according to [14]. The graphene volume frac- 
tion in the proposed RVE model was varied by changing 
the size of the graphene sheet, leading to graphene vol- 
ume fractions ranging from 2% to 6%. Filler volume frac- 
tions of reasonable and practical magnitude were thus 
studied, omitting the agglomeration effects in graphene 
nanocomposites with high filler content. 

The results obtained from the developed RVE structures 
were plotted in Figure 4. These data show a significant 



improvement in the buckling performance of nanocompo- 
sites under compressive loading. The buckling strength of 
neat epoxy was herein considered as the reference level. 
Buckling strengths that are calculated from the direct and 
indirect approach are in good agreement, which validates 
the proposed numerical technique. In the current work, 
up to 26% improvement in the buckling strength of epoxy 
was estimated for only 6% volume fraction of graphene. 

Conclusions 

In this study a representative volume element method 
was successfully employed to investigate the buckling 
phenomenon in graphene/polymer nanocomposites, where 
graphene was assumed to be uniformly distributed. Gra- 
phene was modeled in the atomistic scale and polymer as a 
continuum. A significantly enhanced buckling strength of 
graphene reinforced polymers was observed as compared 
to neat polymer, i.e., buckling strength of graphene/poly- 
mer nanocomposite improved by 26% with only 6% filler 
volume fraction. 
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